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ee RODUCT ON 


A program 1S underway at the Naval Postgraduate School's 
(NPS) Turbopropulsion Laboratory to evaluate the wave rotor 
concept. The wave rotor, operating aS a component ina gas 
turbine engine, uses unsteady wave propagation in tube-like 
passages to compress incoming air before it goes to the com- 
bustion chamber. The combustion chamber output is routed back 
to the rotor to create the unsteady waves. Thus the rather 
Pample rotor, with “partial admission” inlet and outlet ports, 
acts both as a compressor and as a turbine. The alternation 
of hot and cool gases through the same hardware aids cooling 
and allows higher operating (cycle) temperatures to be used. 
The interaction of the hot and cool gases within the passages 
of the wave rotor through the wave patterns that are created 
pose a difficult problem. The work underway aims to develop 
and validate preliminary design and performance analysis tools. 

One of the tools neededis a computer program which can 
be used to construct wave interactions in a design process 
and then accurately predict the performance of the design. 
Met! ©fficient and accurate methods of solving the full 
Navier-Stokes equations for unsteady turbulent flows are 
developed, the design program will be based initially on the 
solution of the unsteady Euler equations and a method devised 


to represent losses. Once the program is developed, it must 


ies 


be verified against experiment. This is the overall goal of 

the present program in which a wave rotor apparatus has been 

assembled and methods of meaSuring the unsteady pressures and 
temperatures are being developed concurrently with the compu- 
tatLonal errore. 

Three different approaches to the solution of the unsteady 
Euler equations were examined in the overall program. First, 
Eidelman developed a two-dimensional code based on the 
Godunov method of solution [Ref. 1] and applied the code 
to examine unsteady wave propagation in ducts [Ref. 2] and 
the process of port opening to wave rotor passages [Ref. 3]. A 
Summary of Eidelman's work is given in Reference 4. While 
the method 1s conservative and does not require the introduc- 
tion of artificial viscosity, the extension from one to two 
dimensions is not rigorous when shock waves are present, and 
computational times with the present code are quite long. 

The extension to include viscous effects would require a separate 
treatment of the boundary layer. 

Second, Mathur developed a ane-dimensional code based on 
the Random Choice Method (RCM) of solution [Ref. 5]. Somewhat 
Similar to the Godunov method, in that the Solution is base 
on solving the Riemann problem within each grid cell at Sach 
time step, the RCM approach results in very sharp discontinui- 
ties which can be tracked easily. This is particularly useful 
in constructing wave rotor cycles, in which the posttlomeas 


the gas-gas interface is equally as important as the position 





of the compression and expansion waves. The code is there- 
fore valuable in the preliminary design process, to examine 
Suitable port arrangements, and the gas properties at the 
ports, for a given task. Unfortunately, an extension to two- 
dimensional and/or viscous flow can not be made rigorously. 

The third approach was followed in the present work. The 
QAZID method for compressible inviscid flow computations 
developed by Verhoff and O'Neil [Ref. 6] was implemented to 
generate a one-dimensional unsteady Euler code with the goal 
of evaluating the suitability of the approach for wave rotor 
applications. Some advantages of the QAZ1D method were recog- 
nized as being the following: 


men The methoaw@is based on the use of characteristics. 
Such methods can model wave propagation accurately. 


2. The use of a natural streamline coordinate system 
eases the difficult task of computing with two and 
three dimensional grids. 


3. The equations are written ina form which allows a 
Straightforward extension to viscous flows. 


4. Codes for computing internal, steady, two-dimensional 
flows in the presence of shocks, and simple internal 
viscous flows, have been generated quickly without 
Significant development problems [Ref. 6]. 

In the work reported herein, a one-dimensional FORTRAN 
code based on the QAZ1D method was developed using the NPS 
IBM370-3033 computer and subsequently exercised on the snock 
tube test problem. In reporting the work, first, in Section 


Il and Appendix A, a complete account is given of the deri- 


vation and non-dimensionalization of the governing equations. 
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In Section III, a FORTRAN code, EULER-1l, is described. Its 
operation on the NPS computer, and thelisting of the code, 

are given in Appendix B and Appendix C, respectively. 

Results of applying the code to the shock tube test problem 
are given in Section IV. Difficulties encountered in the 
implementation of the method and additional comments are given 


in Section V, and concluSions are given in Section VI. 
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A. OVERVIEW 

The QAZ1D method uses Riemann-like variables with a modi- 
fied entropy term in their definitions, to express the Euler 
equations in a natural streamline coordinate system. The 
resulting partial differential equations (PDE), when cast 
along the characteristic trajectories in the space-time 
domain, reduce to a system of ordinary differential equations 
(ODE) which may be solved explicitly. The advantage in uSing 
the modified Riemann variables is that they are less affected 
by discontinuities in the flow than are the standard Riemann 
variables. Since the equations are not valid across discon- 
tinuities which cause irreversible losses, such discontinui- 
ties must be located and treated with special logic in the 
numerical formulation. 

In the following presentation, the reader's familiarity 


Paeema the method of characteristics is assumed. 


By. DEVELOPMENT OF THE EQUATIONS 
Pils Section presents the governing equations and outlines 
their development. A rigorous derivation of the equations 
is presented in Appendix A. 
1. The Coordinate System 
The natural streamline coordinate system (s,n,m), 15 


iene memege A-leicclative to a fixed rectangular cartesian 


Ly 


system (x,y,z). The (s,n,m) system is a right-hand ocrthegene 
system which undergoes curvilinear translation as it moves 
with a fluid particle along a streamline. The coordinate 
system is described in detail in Appendix A. 
2. Variables 
The modified Riemann variables, or "extended Riemann 


variables" [Ref. 6:p. 1], are defined as 


q + AS 


Le) 
tl 


ee SNS, 


where g is the velocity magnitude, A is the speed of sound 
and S is the modified entropy defined in terms of pressure 


and density as 


R 
G , ‘ 
= oo - ] ) 
S 7 (ya) Uae In(P/p )] (2) 
The modified entropy relation is the result of defining the 


entropy change, dS, to be given by 
AG. 20 ae ee (3) 
where dQ, is the heat required to be added in a reversible 


process between the same end states, T is the temperature and 


y is the ratio of specific heats. 
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Stream tube of variable cross sectional area, 


coordinates, 


tives, 


ae) 
at 


where oa 


a 


wise spatial dimension, 


Conservation of Mass 
By applying the continuity equation to a differential 
in natural. 


and ignoring all third order and higher deriva- 


one obtains 


aq 


as 


d0 


| 98 3 
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ee COs OM 0 


rote velocity, 


O 1s the density, s is the stream- 


and 8 and 6 are the flow angles as 


defined in Fig. A-l of Appendix A. 


4. 


law in natural coordinates, 


Conservation of Momentum 
Applying the vector form of the momentum conservation 


the equation of motion for inviscid 


Flow becomes 


@ueside 


changes 


are permitted), 
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The Conservation of Energy 

In the absence of friction and heat conduction, and 
of discontinuities (through which irreversible 
consistent with mass, momentum and energy conservation 


energy conservation is equivalent to the 
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statement that the entropy of a fluid particle does not change. 


In the natural coordinate system, therefore 


dS _ 328 aS 
ae at as 


Equation (6) will not be valid across shock waves or contact 
surfaces between gases having different states. 
6. Transformation to a USeful Form 
Equations (4) and (5) are first expressed in terms 
of the primary variables q, A, S, and the flow angles. Using 


the definition of sound speed in a perfect gas 














2 
A = ¥EVe (7) 
Eq. (4) becomes 
oA 9A , Yol , 3g yr1 , 98 9b, LX 
ae Ol ig Oe ge 1S 8 iO er ee m2) 
With Eq. (3), the equation of state for a perfect gas, and 
the first law of thermodynamics, Eq. (5) gives 
ag oq 2A dA 2 oe 
fe) ee) | 2 ae ee Se 
cL poo, Jee AC cone (9) 
oe acs yq an 
2 
CDi q ae A ceeey Je) 
at aS VAG COS 0) cm 
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Equations (6), (8) and (9) constitute the system of P.D.E.'s 
which describe the unsteady isentropic flow of a perfect gas 
in natural coordinates. Since third order terms were 
neglected in the conservation of mass, the system is only 
accurate to second order. 

The system of P.D.E.'S 1S now transformed into a system 
of O.D.E.'s along the characteristic trajectories. In general, 


may 1S a function of (s,t), then 


OW Ow 


aw = ine ds AP FE at 
or 
dw ow ,dS aw 
SO ddl cena ieee 1 
at joc! Mae (Le) 


where (ds/dt) describes the "characteristic" direction in the 
(s,t) plane along which Eq. (10) gives the rate of change of 


the parameter w. After the appropriate algebra and the intro- 


duction of Eq. (1), the final system of equations becomes 
3Q DOF = it 2 a 
aa (qtA) xe = 5 A(S yal! bas ea A)] 
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The characteristic directions in the space-time plane are 


clearly at+A, q-A, and q. 


C.. SOLOUTPON METHOS 


Equation (ll) may be expressed as 


OW OW _ 
ot 


where the vectors are defined as 
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PeeetajeCeOry in the space-time domain 1s illustrated in Fig. l, 
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Figure 1. Solution Procedure in the Space-Time Domain 


which shows an infinitesimal interval of space between two 
Ppmedes' of the computational mesh. 
For the change in w from A to B along the trajectory 


mweeem Slope A, 
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where 6w denotes the change due to time at a fixed location 
and Aw denotes the change due to displacement at a fixed 
time, for the characteristic trajectory. The essence of the 
solution procedure is to calculate the change in the varia- 
bles at each spatial node during the differential time 
interval, so that the step to the next time level can be 


made. In other words, 6w must be calculated at each node 


using 


7 = er ot oe 
6w = -Aw + il 


N 
Oy 
ct 
= 
& 


Since all the information at time level t is known, the 
position of A can be determined by an iterative procedure 
based on \} and Aw can be calculated by interpolation. —f1— 


line integral can be evaluated by transforming the integral 


to a purely spatial integral using \ = ds/fdt. Thus 
t+5t Z B 7 2 Zz 

f Ze Oe = if 7 ds = X ds 
fe A s,tAs 


and the integral can be evaluated by any one of several 
numerical methods. Thus, the variables at each node can be 


updated in time using Eq. (14) and the process repeated. 
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See OLSCONTINUITIES IN THE FLOW 

There are several types of discontinuities that must be 
considered. They are gradient discontinuities, contact 
discontinuities and shock waves. 

Gradient discontinuities are characteristic of the head 
and tail of rarefraction waves, the collision of two shocks and 
the interaction of a shock and a contact discontinuity. 

Across such a gradient, the derivatives of velocity, pressure 
and sound speed are discontinuous. 

A contact discontinuity is caused by the interaction of 
two shocks of opposite family and when a shock overtakes a 
shock of the same family. Entropy and sound speed are discon- 
eamueus at a contact discontinuity. 

Across shock waves, velocity, pressure, density and 
entropy are discontinuous [Ref. 7]. 

pence Eq. (1) 1s not valid across discontinuities, addi- 
tional logic is required in the numerical procedure to model 
flows which contain discontinuities. The method presented in 
Reference 6 for making a correction in the case of shock waves 
will give good accuracy in problems where the solution is 
Converging to a steady state condition but will not give 
accurate results during the transient portion of such problems 
or for problems with only unsteady solutions. In the case of 
applications to the wave rotor, it would certainly be necessary 


to Know the locations of the contact surfaces between the 


hot and cool gases as well as the locations of the shock and 
expansion waves. 

The methods used in the present work to correct for dis- 
continuities are those of Moretti [Ref. 7]. The present 
discussion will be limited to the treatment of shock waves. 
The method makes use of the analytical relationship between 
the change in the Riemann variables across a shock and the 
incoming Mach number relative to the shock wave (W) which is 
illustrated in Figs. 2 and 3. The relation 1S used to 
determine the shock speed and to transform the problem to 
a steady case which can be handled using normal shock rela- 
tions. For the situation depicted am fig. 92 GF fansioern 
propagating to the right with velocity ie IntO alr wisem 
Ve lLoGieEy dp: and with the high pressure side to the left, 


where A denotes the left side of the shock and B the right, 


u = (qa — a 
WwW = = = SS = (15) 
B B 


If the pressure and density are non-dimensionalized by the 
values on the low pressure side of the shock, the change in 


the extended Riemann variable Q is given by 


= 2 
a 2B 2(We-1) , 2 "a4, 
An (y+1)W f= A 
A 2) y 
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- bs }2£nt [— we - (——) ] (————] ! (16) 
AS iy = 1s) (+) y+1 (y +1) W2 


26 


QA-QB/AB 


2.0 


ol qn a poi a B 


unsteady B steady B 


Figure 2. Shock Wave with High Pressure to the Left 
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Figure 3. Shock Wave with High Pressure to the Right 
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where 
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This exact relationship can be approximated by the polynomial 


a B = -2.7574 + 3.1573W - 0.2863W~ (17) 





as illustrated in Fig. 4 over a shock strength range of 1.0 
to 4.0. It should be noted that if the high pressure side 


were to the right, as in Fig. 3, Eq. (15) would become 


and the left side of Ea. (16) would become 


The procedure then, is to measure the change in Q across an 
interval where a shock is known to exist. This value, call 
ie AQ”! 1s used in Eq. (17) to get an approximation for the 
corresponding value of W and then the exact value of AQ, call 
ike AQer is calculated using Eq. (16). If the exact value of 
AQ 1S not equal to the value measured across the shock, a 


new value, “40, 1s calculated vaccora inate 
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AQ. / = AQ; + (AQ, 7 49.) 
and entered into Eq. (17) to obtain another value of W. When 
mien exact Value of AO calculated from Eq. (16) equals the 

measured value of AQ, W is known to be correct and the normal 
shock relations can be used to calculate the values at node A. 
Also, ee can be calculated from Eq. (15) and used to track 


the shock during the next time interval. 


Zo 


TII. DESCREPTIONSOES FORTRAN PEL OGRAM UGE ae 


A. MACHINE AND LANGUAGE 

The EULER-1 program is written in VS FORTRAN and runs on 
an IBM 3033, System 370 computer, however the code is simple 
and small enough to enter and run on a mini or micro computer 
in a Basic language if one is willing to accept significantly 
longer run times. Table 1 at the end of this section 1s a 
summary of the editable parameters and their effect on the 


PrEocram- 


Bx CONVENTIONS AND “BASTC  SLRUCTURE 

EULER-1l is a first-order one-dimensional code using an 
evenly spaced numerical grid. All values are double precision 
except those used in the graphics routines which are rounded 
tO. sing lev pree ison. 

Each subroutine has its own variables space. In other 
words, variables are not shared in common throughout the 
program but must be passed to the subroutine being called 
by the calling routine. Tn all cases, however, the variables 
have the same name in both the called and the calling routine 
so there should be no confusion. This was done so that 
arrays could be dimensioned at execution time. 

The program, depicted in Fig: 5, 1S Structured sarounaa 
main routine which serves as a user input area, sets up the 


problem and calls five subroutines to solve the problem and 
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output the solution. The five subroutines are mele ca RAK , * 
"SWEEP," "JUMP" and then one of several output routines depend- 
ing on user desires. Output options include two types of 

Seep nical displays, “PLOT" and "EXACT," and one tabular 

hesting routine, “LIST.” When "PLOT"”™ is Selected, "BORDER" 
Bemevomatically called to set up the plotting area. There are 
two other subroutines, "RSHOCK" and "LSHOCK," which are called 


by "SWEEP" as needed. 


fe oUBROUTINE DESCRIPTIONS 

In general, each subroutine begins with a heading followed 
by variable definitions where appropriate, variable declara- 
tions and array dimensioning. Most variables are defined in 
the "MAIN" routine and only those variables which were not 
are defined in the subroutines where they are used. 

1. The "MAIN" Routine 

This routine forms the main structure of the program 
and includes a heading, an extensive list of variable defini- 
tions, a user input area, the necessary statements to make 
the initial value assignments to variables, and the call 
Statements for the various subroutines. 

The input area is those lines of the program (145-175) 
where the user edits the program to establish the initial 
conditions, mesh size, termination criteria and output options. 

ici nieetda | coOncdicloms whach May be modified are the 
pressure, temperature and density ratios across the diaphragm, 


the initial velocity of the fluid on each side of the 


diaphragm, and the value of gamma (which must be the same for 
both sides). 

All velocities are non-dimensionalized by the initial 
sound speed on the low pressure side of the diaphragm and 
pressures and densities are non-dimensionalized by their 
initial values on the low pressure side of the diaphragm. 

The mesh size may be set at any odd number and the 
arrays in the user input area must be dimensioned as such. 

The criteria for program termination is the number of 
time steps computed, which the user selects. 

The output options are determined by the value of the 
variable GRAPHS. A value of zero results ina call to "LIST" 
which writes to file 9 on the::user's permanent disk, a tabular 
listing of the variable arrays and discontinuity Focarian 
A value of 1 results in acall to "BORDER" and "“PLOD” “white 
create a plot of the pressure, density, velocity and entropy 
distributions as in Fig. 7. A value of 2 results ina call 
to "EXACT" which creates a plot of the density distribution 
compared to the exact solution as in Fig. 8. The exacuevama 
to be plotted must be entered in the user input area, they are 
not calculated by the program. A more detailed description 
of the various outputs is found under the appropriate subrou- 
tine description. The frequency with which output is created 
is controlled by the variable SKIP which is the number of time 


steps between calls to output routines. 


Zoe thie RIMES Routine 
The maximum allowable time step iS governed by the 
CFL condition. Simply stated, this means that the time step 
must be small enough so that the characteristic trajectories 
remain within one spatial interval during the time interval. 


The minimum time step is calculated by computing 


Delt = H/ABS[Q+A] 


at every node, where H is the spatial interval, and selecting 
the minimum value of Delt. 
Ee fhe STRAK" Routine 

Shock locations at the unknown time level are deter- 
mined by computing the shock speed at the known time level, 
as outlined in Section 2, and multiplying by the time inter- 
val computed by "TIME." The time step 1S reduced, if neces- 
sary, to limit shock travel to one spatial interval. The node 
immediately to the right of the shock (upstream) is flagged 
fmemelater use by "SWEEP" and "JUMP." 

In calculating the shock speed, the assumption is made 
that the conditions immediately adjacent to the shock are the 
Same as the conditions at the nodes to the left and right of 
the shock. 

4. The "SWEEP" Routine 
The "SWEEP" routine makes the necessary calculations 


Bemoolve EQ-s(4). This involves interpolation at the known 


SiS, 


time level for the values of Aw, Calculatdentol gland 
integration of Zz for each mmoder 

To calculate Aw, the assumption is made that tue 
characteristic trajectories are straight lines. An initial 
guess of the characteristic slope (A) is made and by forcing 
the characteristic to pass through the node at point B in 
Fig. 1, As = X} xdelt. The values of q and A at point A are 
found by interpolation and the slope of the characteristic 
is then calculated using the values of gq and A at point A. 
This slope is compared with the initial guess and if the two 
are not in agreement, the procedure is iterated until thev 
converge. Once As is accurately determined, the values of Aw | 
can be calculated by interpolation. Linear interpolation is 
used in the present version of EULER-l. This procedure is 


carried out for each Characterisete at ecachmmeqce, 


Gam ene 
LL 
S 
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A deviation from the above procedure is necessary in 
the case where a shock exists in the spatial interval to the 


left or right of the current node. When a shock exists to 
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e@enrighnt, and the flow 1s subsonic, aS in Fig. 6, care must 
be taken not to interpolate for point 3 based on the slope 
between points I+l and IT, which would not be accurate due 
to the discontinuity. In such a case "RSHOCK" is called and 
the interpolation is based on the slope between point T and 
I-l. Similarly, when the shock is to the left of the current 
node, "LSHOCK" is called and the interpolation is based on 
the slope between I+1l and I. 

tre calculation of z includes the calculation of 
Spatial derivatives of q and A for characteristics 1 and 3 
of Fig. 6. In general, the derivatives associated with 
characteristic 3 are forward differenced and those associated 
with characteristic 1 are backward differenced in keeping 
with the principle of domain of dependence. If the flow is 
Supersonic or if a shock exists to the right of the current 
node, all derivatives are backward differenced. If a shock 


Sxists to the left of the current node, all derivatives are 


forward differenced. Once the derivatives are known, z is 
calculated from Eq. (11) using the average value of A between 
Seances | and tf or 3 and I aS appropriate. Since the deriva- 


Meves are linear, this results in an average value of Zz 
over the same interval. 

The integration of z is transformed from a time to a 
Spatial integration as described in Section II and the 
trapezoidal rule is used to carry out the integration. In 
BFULER-1, this has been done in one step using the average 
value of z described above. 


Su 


Equation (14) is then solved by addition of Aw ane 
the results of the integration, and the resulting 6w is 
stored until all nodes are calculated. After all nodes have 
been calculated, the variable arrays (w) are updated for the 
next time interval. 
5. The "LSHOCK" Routine 
As discussed above, the "LSHOCK" routine modifies 
the basic EULER-1 procedure at an interior node when a shock 
exists to the left of the node. Interpolation of guantities 
to the left of the node are computed based on the derivatives 
of the quantities to the right of the node. Although this 
assumes that the derivatives do not change between adjacent 
spatial intervals, it is necessary to avoid taking derivatives 
across discontinuities in the flow. 
6. The “RSHOCKY Routine 
Similar to "LSHOCK," "RSHOCK" bases interpolations 
of quantities to the right of the node on the derivatives of 
the quantities to the left of the node when a shock exists 
tO the Fright Of Ehe node. 
ie The “JUMP” Routine 
The "JUMP" routine is used to calculate the conditions 
downstream of a stationary shock as described in Section If. 
If a shock is known to have crossed a spatial node during a 
time interval, which is known once "TRAK”" has been cal ledmaen 
that time interval, the entire "SWEEP" sequence is skipped for 


the node which was crossed by the shock and the conditions at 


38 


that node are determined using the normal shock relations. Note 
that the node in question is the node downstream of the sta- 
tionary shock. As in "TRACK," it 1s assumed that the condi- 
tions at the nodes upstream and downstream of the shock are 
the same as the conditions immediately adjacent to the shock. 
8. The Output Routines 

There are four subroutines which produce various types 
of output as the user desires. "BORDER" and "PLOT" produce 
a graphical presentation of the pressure, density, velocity 
and entropy distributions at selected time levels as seen 
for example, in Fig. 7. "EXACT" produces a plot of the den- 
Sity distribution at one selected time interval and compares 
it with the exact solution as shown in Fig. 8. At selected 
time levels, "LIST" produces a tabular listing of the Riemann 
variables, modified entropy, pressure, density and velocity 
distributions, elapsed time, shock speed and discontinuity 
locations. The listing is written to the user's permanent 
storage disk. 

When the value of GRAPHS is set equal to one in the 
"MAIN" routine, "BORDER" is called once to set up the plot 
axis, labels, and headings. "PLOT" is called every SKIP 
time steps to draw the four distribution curves. 

When the value of GRAPHS is set equal to two, "EXACT" 
is called every SKIP time steps to plot the density distri- 
bution compared with the exact solution computed at six 


points of interest. The points are the two end points, which 


ary, 


are simply the initial conditions, the head and tail omen 
rarefraction wave, the point just left of the contact surface 
and the point just left of the shock. The spatial locations 
of these points are computed by "EXACT" based on the elapsed 
time and the known values of the wave velocities entered in 
the "MAIN" routine. The exact values of the densities at these 
points must also be entered in the "MAIN" routine. 

When the value of GRAPHS is set equal to zero, the 
tabular listing as described above is sent to the user's 


disk. No graphical output is created. 


D. CAPABILITIES AND LIMITATIONS 

The present version of EULER-1 is set up to solve a shock 
tube problem with a single centered diaphragm. Boundary 
conditions for the ends of the tube have not been incorporated 
so the problem must be stopped before the waves reach the end 
of the tube. 

The left side of the diaphragm is the high pressure Side. 
The program will] not run with the right side as the high 
pressure side without changes to some of the shock corpeemiag 
[Weye alex 

The user may select any odd number of grid points Piumiees 
only by the amount of memory available. 

The program tracks shock waves and makes shock jump calcu- 
lations at the appropriate locations but the program can also 
run without the shock tracking feature and jump calculations 
if so desired with some smearing of the shock discontinuity 
and complete loss of entropy change across the shock. 
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TABLE Ll 


LIST OF EDITABLE PARAMETERS 


Parameter Punee Lon 
N Numer OfGmaaid POInts 
M Controls which discontinuities are 

tracked 
Me oaocksS Omly 
2 = Shocks and contact surfaces 

GRAPHS Coueeeils the foOrmof the output 
UV —-aemlar- Listing 
l = Pressure, density, velocity, 

entropy plot 
2 = Exact solution comparison for 
density 

SKIP Number of time steps between output 
calls 

iS LOP Number of time steps calculated 

eke 1 Initial temperature ratio 

PRI Initial pressure ratio 

DRI Initial density ratio 

Oi Tied veloc hey wlert of thevdlaphragm 

ORI Initial velocity right of the diaphragm 

G Ratio of specific heats 


4l 


LV eS teers 


EULER-1l was tested on the Riemann shock tube problem. 

The initial conditions and the courseness of the computa- 
tional mesh were varied. Additionally, one run was made 
without the shock tracking and correction features. 

Figures 7 and 8 show the results for initial pressure 
and density ratios of 5, uniform temperature and with the air 
initially at rest, uSing a grid of 101 points. The exace 
analytically predicted conditions are also shown in Fig. 8. 
It can be seen that all wave velocities are correctly computed 
and the shock is defined within one interval. The rarefrac- 
tion waves are Slightly smeared and the contact surface is 
greatly smeared. A slight transient instability to the right 
of the diaphragm is noticeable, particularly in the plots 
of pressure and velocity. 

Figures 9 and 10 show results for initial pressure and 
density ratios of 1.3. The shock is still well defined in 
the correct interval and the contact surface is not as badly 
smeared aS in Fig. 8. There is a loss of accuracy however, 
with respect to the tail of the rarefraction wave. 

Figure 1l shows the results of the original problem with 
a grid of 51 points. AS might be expected with a coarser 
mesh, the discontinuities are less sharply defined although 


it is clear that the wave velocities are accurately computed. 


Figure 12 shows the results for the Spiga yproblem with 
a grid of 901 points. It can be seen that the BULER-1 
solution closely approximates the exact solution. 

cures 13 and 14 show the result of running the original 
problem without the shock tracking and correcting features. 
Note that the shock position is incorrectly computed and it 
1s smeared slightly. Also note in Fig. 13, the lack of any 
change in modified entropy across the shock. 

Run time for the 101 point mesh is 0.00049 seconds per 
node per time step. The results for Fig. 7 took 50 time 
steps for a total time of 2.46 seconds. For the 901 point 
fesm che run time decreased to 0.00023 seconds per node per 
time step due to the less frequent use of the shock tracking 
and correction routines per node calculated. The results for 


meg t2 took 470 time steps for a total time of 71.85 seconds. 
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A. RESULTS 

The transient instability noted in the results was not 
investigated since the amplitude was small and the effect 
would be of no consequence in most applications. 

At low pressure ratios, the velocity of propagation of 
the tail of the rarefraction wave approached sound speed 
and, in these cases, the EULER-1 code did not accurately 
compute its location. The reason for this was not investi- 
gated here; however, due to the Jow pressure ratios at which 
wave rotors can operate, this is of interest and should be 


investigated in future work. 


Ba SPECIAL CONSTRPERATIONS 
1. Modified Entropy 

Tn the QAZ1D method, and subsequently, in the EULER-1L 
code, the variable S, which here has been referred to as the 
“modified entropy" and which is referred to in Reference 6 as 
Simply the "entropy," does not behave thermodynamically as 
one would expect entropy to behave. As a fluid crosses a 
shock wave, the modified entropy of the fluid decreases. 
This is the result of the definition of Eq. (3), repeated here, 


that 


SS) == = 


a2 


Miesuse Of the word “entropy” and symbol "S" for this varia- 
ble 1s confusing since the sign of the change in "S" is in 
direct conflict with well-established conventions in thermo- 
dynamics. For example, Corollary 6 to the qe Law of 
Thermodynamics [Ref. 8:p. 88] would have to be changed to 
read, "The ‘modified entropy' of an isolated system decreases 
Or in the limit remains constant." At a minimum perhaps, the 
eymool S to denote the “modified (negative and scaled) entropy" 
would be helpful. In practice of course, the change in the 
conventionally defined entropy can be recovered from the change 
in the modified entropy, Simply by multiplying by (-y). 

2. Moretti's Methods 

Incorporating Moretti's methods of handling discon- 
tinuities [Ref. 7] into the EULER-1l code required special 
care. 

There is a fundamental difference in procedure in 
that the EULER-1 code is based on the high pressure side 
being on the left whereas Moretti's formulations are all 
based on the high pressure side being on the right. Although 
essentially a matter of bookkeeping, it is an area of poten- 
fal confusion. 

Another difference is that the Riemann variables used 
by Morreti are not the Riemann variables used in the QAZ1D 
Pmeenoed: This Can lead to difficulty. ime tact. 1n the 
omtedt lOnmeoreMorectti's Shock tracking scheme to the EULER-1 


code, the pressures and densities have to be non-dimensionalized 


ay 


by their values on the low pressure side of the shock in order 


to obtain Eq. (14), and for the procedure to work. 


C. SUITABILITY FOR WAVER O1LORe=e sh terrrerls 

Further work is necessary before the suitability of the 
QAZ1D method for wave rotor applications can be fully evalu- 
ated. The following areas need to be addressed. 

].  BeoundasyeGondicions 

In the EULER-1 code, boundary conditions at the ends 
of the passage have not been incorporated. The code calculates 
the changes of the interior nodes only, skipping the two end 
nodes. In particular, solid wall boundary condie1ons jane 
open-end boundary conditions must be incorporated. No particu- 
lar problems are expected in the boundarv conditions themselves. 
However, additional logic may be necessary in the handling Ot 
discontinuities, such as the reversal of conditions to the 
left and right of a discontinuity atter reflection f£rome 
pounce The additional logic, which may be considerable, 
LS warranted by the simplicity of the QAZ1D method. 

Z. €oOntact Disconeinuie, 

No special attention is given to the contact discon- 
tinuity in the EULER-1]. code and hence the discontinuity 1s 
smeared. In a wave rotor application, this discontinuity 
will have to be sharp. Moretti's methods again appear to he 
applicable here and the additional logic needs only to be 


formulated and incorporated inta the EULER-1 code. 
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3. Quasi-One-Dimensional Modeling 
EULER-l is a one-dimensional code. The solution of 
flows in passages of varying area may be desirable for some 
wave rotor applications. Such problems can be solved ina 
quasi-one-dimensional manner by the addition of the appro- 
priate area variation term to the equations, and the need to 
solve the fully two-dimensional equations is avoided. The 
area variation must be incorporated into the EULER-1 code and 
the code tested on quaSi-one-dimensional problems. 
4. Other Potential Extensions 
Another capability which can be incorporated is the 
ability to handle two gases with different specific heat 
ratios. Eventually, the code should be extended also to a 
second order, one-dimensional version and then to a first 


Order, two-dimensional version. 
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Vi. “@NCLUS#ENS 


The development of the EULER equations in the natural 
streamline coordinates using extended Riemann variables was 
reviewed in detail. The QAZ1D solution method, with the 
addition of shock tracking features, was implemented ina 
first order, one-dimensional FORTRAN code with the intent of 
evaluating the method's suitability for wave rotor applica- 
tions. The code (EULER-1) was tested on the shock tube 
problem with good results. The incorporation of boundary 
conditions, an improvement in the contact discontinuity 
definition and the addition of an area variation term for 
quasi-one-dimensional modeling are considered to be neces- 
sary before the QAZ1D method can be accepted as being suitable 
for wave rotor applications. The additional logic requieee 
for these extensions may be considerable but the development 
1S warranted by the overall simplicity of the QAZ1D method, 
and its Straightforward extension to viscous, multi-dimensional 


Flow modeling. 
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APPENDIX A 
DERIVATION OF THE GOVERNING EBOUATIONS FOR 
THREE-DIMENSIONAL INVISCID FLOW (WITH ARFA CHANGE) 

The governing equations of a 3-D, inviscid, compressible, 
unsteady flow are derived here in a natural streamline 
coordinate system. The equations are then recast along 
characteristic trajectories in the (s,t) plane and expressed 
in extended Riemann variables, reducing the system of par- 
eral differential equations to a system of ordinary differ- 


ential equations. 


A. DEFINITION OF VARIABLES 


A Cross-sectional area of a differential stream 
tube 

A Speed of sound 

Cy Specific heat at constant pressure 

ee Specific heat at constant volume 

da Normal vector for differential area of control 
Surface 

e Specific internal energy 

1 UME ctormsminimiier sie sand im directions 

Sin 

P Static pressure 

QO Modified Riemann variable 

Op Reversible heat transferred 

q Velocity magnitude 

xr Position vector 


a7 


R Modified Riemann variable 


Ra Gas constant 

S Modified entropy 

T Static temperature 

Ic Time 

u Velocity relative to a standing shock wave 
V Specific volume 

vol Differential element of volume 

V Velocity vector 


Ww Incoming Mach number relative to a standing 
shock wave 


Q Density 

Y Ratio of specific heats 

8,0 Flow angles with respect to reference coordinate 
planes . 

V VEGeer (epee gor 


Be THE COORDINATE aS io FEM 

The natural streamline coordinate system (s,n,m) 1S 
Shown with respect to a fixed rectangular cartesian system 
(x,y,Z) in Fig. Al. The system is a curvilinearly transi@eagee 
right hand orthogonal system which translates with a fluid 


particle along a streamline such that the 's" coordinate 15 


measured in the direction of the flow. The 'y' coordinate 


direction always lies in the plane defined by the 's' coor- 
dinate direction and the fixed ‘y' axis. The '"m'’ direction 
normal to the (s,n) coordinate surface. The flow angles, 


and ¢, are defined as shown in Fig. Al. 
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Figure Al. The Natural Coordinate System 


C. VECTOR OPERATORS IN THE NATURAL COORDINATE SYSTEM 


|. a. 





—> 
Licino wetenpOstitem vector Of a fluid particle 


ie natural coordinates, Eq. (Al) becomes 


“nN ”N 


[ids + i dn + idm] -{( ee staal ate + ( on 
ee ne mo ae ig ek Lae 
0s on am 


By inspection, 


2 




















ine tk Lae das one a 
cae 5 Iss ‘i yee | Ayo om (A2) 
2. 
Since 
a — le 
using Eq. (A2) 
~~ re) 
VVC) = q gt (A2) 
ae Vey 
From Eq. (A2) with V = Gre 
V-V = ne (Aa 
4. (V°V)V 
From Eq. (A3) with V = oie 
ar Te JV = S = oq. S 
a is 4-95 — Slsate ee aS Sse) 


The change in the unit vector ne with respect to smi 
derived below and is illustrated in Figs. A2 through A6. (A 
three-dimensional model is very helpful in visualizing the 
vector geometry.) 

Figure A2 illustrates the natural coordinate system 


translating along a streamline from point A to point B in 
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Figure A2. Unit Vectors Moving Along A Streamline 


| 
| a N | 
f_- ou 
e e E 
al Ghee Say 3), [Meera Oosii ton Onn tee vectors from Points 
g perp e 


A & B 


61 





Figure As. X=2eP lone 


Figure Ad. Plane Formed by B Vector and Y¥ Ax1s 





oe \ 


Pigjure Ah. Plane Porméed by Al Veeceor ania = 


W 


[UO 


the inertial cartesian coordinate system. Figure A3 is the 
Superposition of the i unit vector from points A and B of 
Fig. A2. The intent here is to express the change in the 
pet vector i in terms of components along the original. 
fenem)) directions shown at point A in Fig. A3. Vector OA 
in Fig, A3 is the original i direction and vector OB is 


A d1 


1. 2 = oa Points cc and Hoare the; projection of points B 


and A in the (x,z) plane respectively. Point D is the projec- 
mio Of point C onto the line OF. 
Figure A4 is the plane formed by the OB vector and the 


Meets. The length of OB is unity by definition. The length 


of BC is sin(6+d@). 
Piece) = "9 Sin 9 cos’dS + cos 6 sin dé 
which, since d@ is a small angle, is 
Sim terao). = s sin Oetecos 6 dé (A6) 
The length of OC is cos(@+d6). 
Goctorag) = Sees oe eos @e - sin 6 sin’ de 
which, since d8@ is a small angle, is 


cos (8+d4) =e Oc. — scanned ae (A7) 


oes) 


Figure A5 is the (x,z) plane in which the angle 94 is 


defined. The La direction is shown at point C. The length 


of CD is 0C sinddé which, with Eq. (A7) is 
[cos 6 = Sin geeae |sinwed 
which, since dd 1s a small angle, is 
G@S 6 doe=—tsialiwc- ademae 


Dropping the second order term 


—— 


CD = cos @ do (A8) 


The length of OD is OC cos d¢ which, since dod is a small angle, 
Se, 

Figure A6 is the plane formed by the OA vectoreane 
the y axis. GD is the projection of BC and since both ese 
vertical lines, BC and GD have the same length, which has 
been determined in Eq. (A6). The — and i directions are 
shown at point A. The lengths to be determined are AF angmeue 
since these are the components of AB in the s and n directions. 
Ror small angles, OD was sShowm Loupe wequahmno OC and GD =sEa 
so that the angle GOD 1S equal to angle BOC = 6+d0@. Angle GOF 
is therefore dé, the length of OG is unity and the length 
of GF is sin d8, which for small angles is d@, and it is in 


the is direction. The lengthasof Ar foes Of sen un Un ies, 
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by definition and OF is cos d6, therefore AF, the component of 
AB in the i direction is equal to 1-cos dé, which, since 
dé is a small angle, is zero. 

The component of AB in the a direction is GF or 


dé i, (AQ ) 


The remaining component of AB in the i, direction is 


ep or 


~~ 


cos 6 dd La, (A10) 


Therefore, the vector AB can be written as 


~ 





a1 i 5 
S = . e 
ae ds = dé 1, + cos 6 d¢ te 
moo ei a . 
Ste Gis 1, + cos 6 ve ds i, 


fee dividing by ds, finally 





o1 ”™ ™ 
SS . ob Op . 
7S eee COS 0 a te (All) 


Note that Fig. A3 can be considered also to describe 


a change in the unit vector 1. deena location, with 


respect to time. In this case, the vector AB is equal to 
d1 
Ss 


sr Siemoncewati hoe (A9) and Eq. ({Al0) 


oS 





O1 an m 
s _ 263 a¢ 
stn | CO SE eer (A12) 
With Eq. (All), Eq. (A5) becomes 
— = Bee. =. 2eou ee 2 od 7 
(WT = aggistad gg i, +d cose 55 i, aio 


D. CONSERVATION OF MASS 


In the quasi-~one-dimensional differential stream tube, 


Shown in Figure A7, where is 1s always in the direction of 
V, 0,4, and the cross sectional area A are known at these meen 


of the element and p and q are assumed constant on any given 


Gross SeGe lon. 


do ds 
h~ 3s 2 


e) 


ol 


ds 
is 2 
dA ds 
bos) 2) 





Figure A7. Differential Stream Tube of Variable Cross 
SeCtlon 


The statement of continuity is 


The change in the mass The netein & lass om 
Within the’ conesol volume = mass across tne 
with respect to time control surface 
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sanee the volume is constant, 


a 
Ou 
7) 

a 
OD 


Q? 
ct 


For the stream tube, 


the left side becomes 


mass enters and leaves only from the 


Q) 
a) 


fo 


as) 
2 


Q) 
| 


dq ds dA ds 
la-3s g/ 'A-gs Z! 


which on expansion becomes 


JA 3 3 
eed 56 ds - gA a es) — 9 A d 


— ds 
S 


the left hand side, 


On dropping the higher order terms and combining with 


00 IL els 00 0q 
weer in ge " S.55 “Ee se 
Equation (Al14) 


with area change. 


In general, 
am terms of q, 


Giatacesd: . 
With reference to Fig. A8, 
sectional area of the 


Stream Elbe scam be written 





JA - odm adn 
re Ck = [dm + Wen ds|][dn + 


vou S 
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ends so that the right hand side can be written as 





is the form of the continuity equation 


which is required to model a flow as being one-dimensional 


however, 


the change in the 


PA) 


must be expressed 


(Fe 5.) 


Figure A8. 


C 
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After expansion, 


fErms .: BCs, Ate 


which, ror smal 


Bac Ole acre Oe, 


>| 
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_l 
eg 





ross Section of a Diverging Differential 
tream Tube 


dividing by A and dropping the higher order 














becomes 
JA 1 (eee 1 oddm 
oe So —_. ( 
aS dn os = dm ds eR) 
ds = = ds ee any 
dn (28 Sips on 
cos (a7 

angles, becomes 

Gdn. —<. See 

adm Mg i 

aS Gsy 2 am AB (Alea 





Figure A9. (s,n) Plane of Diverging Differential Stream Tube 





Figure AlO. (s,m) Plane of Diverging Differential Stream 
Tube 





pagure All. cn) elastomers sm Py rOoTectlon onto the 
(Ce 27)) Pilane 


Oe) 


With reference to Big 2ecrl AB = CD, and 
CD = OD sin (x2 dm] (Al19) 


OD = ee (A20) 


cos [<2 dm | 
m 


OC “= dsecose ce (AZ is 


Combining Eqs. (Al18) through (A21) and knowing that a dm 


is a small angle, 





= ao 
cos 6 ve (A22) 
Equation (Al6), with Egs. (A117) and (A22), becomes 
allen 7 COSkE a0 (A2 3) 


With Eq. (A23), the general form of the continuity equation 


Im natLluras. @Ceordinates aus 
dg 09 do 
ote, —— — ( 
5 + q oe + O75 + egla— + cos 86 5 ] 0 A2aD 


lade CONSERVATION OF MOMENTUM 
The statement of conservation for an arbitrary control 
volume which is fixed in the reference frame can be written 


as 


70 


The change in the The net influx The net force 


momentum of the fluid} _ | of momentum , |on the con- 
in the control volume across the con- trol surface 
with respect to time trol surface 
The net 
i body 
force 
on the fluid 


Assuming the fluid is inviscid, the only forces acting on the 
control surface are pressure forces. In the absence of any 
body forces (gravity, electromagnetic, etc.), the statement 


becomes 


0 


aE oievedvol! = = {{ Vov-da - ‘ine aa. 0 (A25) 


With Gauss' Theorem, Eq. (A25) becomes 


fff <[V]avol A {f{fv- lVvpVJdvol + {ff VP dvol = 0 
Or 

fff zeov) + V'[Vev] + VP}dvol = 0 (A26) 
Since Eq. (A26) is true for any volume, the integrand must 


be zero, thus 


[pV] + V-[Vovi + VP = 0 (AT) 


Ok 


Expanding the first term of Eq. (A27) and using the vector 
identity 
V-(VoeV] = V{[V'pv] + p(V°V)V 
Eq. (17) becomes 
0 = + vise + V-pv} + O(V°¥)V + YP = 6 (A28) 
The second term of Eq. (A28) is zero by the conservation of 
mass so that Eq. (A28) reduces to 
aV = = 
fe E +O (V = V Veneer e 0 (A29) 
Using Eqs. (A2), (A5), (Al1L) and (Al2), collecting termsemaae 
equating each vector component to zero, Eq. (A29) becomes 
. 3g og , oF 
1.lp ot ye ds ss 7S. (A30) 
- Oke Zea dP 
1, [eq ra + pg ne ate an) 0 (A31) 
‘ 2 dm , OP 
La! Oq cos 0 Tt FelG eo S.9 ae RE eae = 0 (A3Z 
Using the perfect gas relation 
i ees 
ie 


a2 


and the identity 


Qo 
ng 


a ene 
as 


o| 
Y) 
| 
rg 


Eqs. (A30) through (A32) become 





2 
3q ee ee one 
a 6 43s UU Us a 
09 ole Ae eal 1 
32° 49s ~ ~ Yq on ion 
ee _A”___ 3 gnp (A35) 
Jt 19s Yqcos@ om 


Woich are in the desired form. 


fee CONSERVATION OF ENERGY 


The conservation of energy for an arbitrary control volume 


1s 

The rate of increase of The rate at which energy 

energy within the con- _ |is entering the control 

trol volume with respect volume across the 

to time boundary 
The net rate of 
work done on the 

-- : S 

fluid at the EO) 
boundary 


Neglecting gravitational potential, the energy per unit mass 


TS 


We 


q? 

a 
ania) 
where e iS the internal energy per unit mass and g is the 
velocity magnitude. For an inviscid fluid in the absence of 
body forces, the only work done on the fluid is pressure work. 


Thus Eq. (A36) is written mathematically as 


Z 2a “a 
ms {ffle+ 16 avers += -{fle +3 -16V-ad - {f{ PV-da (A37) 


Using Gauss' Theorem, Eq. (A37) iS written as 


2 
(ffate + Zio avol 


” 


2 
- [fJv- tle +> Vlavol 


- {ff V-[PV]dvol 


Oe 
3 q? q* = = 
api le ge osy niel +1 ee ev = Bae 
Expanding, 
20te 42 + 0 Sa + (pV) V+ Le a + [e Ey (9- (0)) 
= ~V-[PV] 
28 
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a 7 : 2 7 2 
le +3] [s2+V-(pV)] + 9 aele +5] 2 (oe) 12 le +35] 


Sy [PV] (A3 8} 


The first term of Eq. (A38) is zero according to the conserva- 


tion of mass. Expanding the remaining terms of Eq. (A38), 
) oe ,) q* = 
Beer eg) ax! 7p) > eV IP) 
or 
D q? 1 — 
pr le pom = - 5 V- [PV] (A39) 


Equation (A39) can be written as 


De Dd ee eg ~ Pog. 
DE P ee x! ou /P) a V) (A40) 


Aw 


With V = qi. in the natural coordinate system, Eq- (A29) 


becomes 
aV A 1 
— ——s = = Pp 
5 a (2 De 5 VI 
or 
DV ] 
—_ = a 1 
DE 5 VP (A 4a) 


PhS 


Taking the dot product of V with Eq. (A411), 


= 
| 
tl 


=(V-9P] (A42) 


and substituting Eq. (A42) into Eq. (A40) 








2 _ 
De D 388) ae DNs Pp _ 
— + —f——|] = —_ - —[( Ve 
De * pela? Te ae (A43) 
Bue 
pv A.) SAEs. epee 
Dt Dt Dts 7 Det 
pa di d1 
pe as 9 ce em 
and using Eq. (All) and Eq. (Al2), 
Bi bars chee . 
ve ee i. + als 1 + cos @ re 1 
7 goes i+ cos 9 22 i J 
S aS m 
Therefore 
= OY a ek Dicien) icles ao 4, 30 
De dis (Be * va os fe) cies ee gc (c1OKS teh les aS! ae 
Z 
Pe so 
ew Das Dis, <2 
Therefore Ea. (A43) reduces to 
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P — 
a 7 = —(V-V) (A44) 
Emaem the conservation of mass 


evan — | -0 


Q>| Ww 
chine 


Or 


which, in the natural coordinate system gives, 

sh + q 22 = -piv-W) 

OM 

oa = (Vy-¥) (A45) 
em@ostituting Eq. (10) into Ea. (9), 


Bes be SG (A46) 


Va 


and the first law of thermodynamics for an elemental mass 


requires that 


_ i) 
dQ, = de + Pd = 
so that 
ul it 
as = - —~Ide + Pd(-) |] 
ail (=) 
= yp ide + 2 do] 


and Eq. (A46) implies that 


mi (9 (A47) 





which states that modified entropv is conserved along streamlines: 


G. MODIFIED ENTROPY EVALUATION 


With the modified entropy defined as 


dQ 
= 
clue = a 
Ore 
B dQ 
ee Pes 
S Gay p= Y -! ; (A48) 


From the first law of thermodynamics 
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dQ = de+P dv (A49) 


where e is the internal energy and v is the specific volume. 


S@eceltuting Eq. (A49) into Eq. (A48), 











al: de +P dv iE de P dv 
e-S, = — f{ =—[ f —=+ f | (A50) 
pe Yaa e a aR z 
For a perfect gas 
Pv = RT and de = Soe ay 
for which Eq. (A500) becomes 
; ; : . B C., olde . B Re af 
A B Yon T A V 
which after integration yields 
1 TR Ve 
S = S, + -[C., 2n — + R &n —] (A51) 
A B Yay Te Vy 
pubpstituting 
= ( _ 
Ro = OY i) 


ieee Bo. (Ad5l) 
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At: V V 
i B B B 
S S =[C_(’n == -&n —) + Cy y Cn —] 
A B Bs A's qT) Va Va 
R R fe 
alt G B BA A VR A 
S S =| ee ) + In —] 
A B yy-l I R, B yr Op 
and since Ra = Ra 
B 
P le 
R 5B B 
S S [2n —— - y Cn —] 
B = 
A Va Pa On 
OG 
S Pp 
oa B 1 Bo PA 
am ie Bey a 
G CG bas PRY Pry 
Or 
S S P Pp 
= _* aay lin S = Qn Py (A52) 
C C VEN aa Pet eeu 
lt the reference condition @2Secnesen tobe 
ae... be l Pa 
Ri =T aCe. ome 
é bs : : A 
EG: becomes 
S P 
B 2 il B 
ee, = = gn — {(A5s) 
Ro ea viv ORY 
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Pee eee tHeER TRANSFORMATIONS 

Equations (A24), (A33) through (A35) and (A47) are in the 
required form. However, some of the variables need to be 
transformed, namely, derivatives of %n P need to be exvressed 
in terms of A and S and derivatives of p must be expressed in 
terms of gq and A. 


Using the first law of thermodynamics expressed as 





qQ = dh - — dP 
p 
Eq. (A48) may be written as 
I 
ol ou he | a VP 
or for the oe component 
a1 ak 
ee eee Pe One (A54) 
' 08 T os ol, os 


Assuming Cy to be constant and using the equation of state 


for a perfect gas Eq. (A54) becomes 





Yas ToS Shee as 
G 
and since ea preter) 
_, 9S = ike (ae a 2 a kn P 
ie (y¥-1) P 8s Rap G as 
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Since Ris a constant it may be moved in or out of the 


derivatives giving 


ie Ooms 
oS 





aa 


a P 
Sp 


7 
tr 
Ilo 


and since y is also constant, similarly 





~ » 22S a ee ee 
Js Ra vy-l1 Py ds JS 
Z 
Using A = yP/o 
2 
as, eo a0 gel 2a eee 
dS Re y-l Ae dS an 
2Ay SA _ 0 %n P 
", oe 25°. eo 
Ge garter az 
and, on rearranging 
0 &’£n BP _ 2y 1 dA Ora 
men ys) bee aaa 


Substituting Eq. (A55) into. Bat Aes | tyres 
S = 
= lm! = 0 (A56) 


and the derivative of @n P has heen removed. 


aid 


Since 


[aye sis 
2 


logarithmic differentiation yields 


For a pure substance 


P = P(0o,S) 


and an isentropic process 





Ps er 
Differentiating 
oP -— GPa YP... Z 
(0S do = 
therefore 
Gee 2. do 
le : 0 


em@lie@seltuting Eq. (A58) into Eq. (A57) gives 


dA Ae ec 2 
2 — = peal aaa = =a 
u £ O u O 


or 
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(A57) 


(A58) 





do = he (ASoao 


Along a streamline 


QU? 
Q> 
OD 


(A60) 


Ou 
mm! 

il 
a 
(er Be: 

an 

re 


Q? 
in 


and since the fluid is the same from streamline to stream 


Eqs. (A59) and Eq. (A60) cam be combined giving 


3 
.) 


= 2p 0A, 5 aA, (A61) 


00 dn 2 eee 
A (y-L)A ot as 


op . 20 
eo dS yo. 


ey 





tT 


Substituting Eq. (A61) into Eq. (A2Z4) yields 


3A oA 


cist ae a) 
as z Z A S ce vA aa 





ao, _ 
+ cos @ am! 0 (A62) 


and the derivatives of 0 have been removed. 


The governing equations are now summarized: 





on Cate hae cy peck imran (eee 20) 2 

; + =e Ax +aA 5 Lor + COS. 6 ae = 0 (A63) 
9g aq , _2A 9A 2 ae 

St isa jee" veep. aeie = * Ce) 
a6 |. 38 _ _ AS 8 tn P (AG65) 
ot =20o. | > vq on 
ek) od AS a <n P 

Sc es (ees 

wea a ae Yq cos@ 4m (A66) 
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aE SS 


Note that Eqs. (A65) through (A67) are in the form 


aX A». 
cee ae 
mewx = X(s,t), then 
ne Ps aX 
dx = 3¢ dt + a5 ds 
and 
Gee — «0% ds 0X 
de. se ae ss UD ey 


mmeme (Aas/dt) 1s the “characteristic" direction in the (s,t) 
plane along which Eq. (A68) holds. Along thisdirection, the 
equations may be transformed from partial to ordinary differ- 
ential equations. What is desirable is to find transformed 
variables in terms of which Eqs. (A63) and (A64) will be in 
the same characteristic form. For this reason the modified 


Riemann variables 


(A69) 


> 
il 


qt AS/Ro 


ve) 
lI 
Q 


= AS/R. (A70) 


a 
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are defined, where S/R. is defined in Eq. (A53). The modified 
Riemann variables can be introduced by manipulation of the 


equations as follows. First, multiplying Eq. (A6G3)) by we 


gives 
OA oA y-1 dq tee 0 > 7 
S5¢ + 9S55 + AS ae 7 GCS a) COS Fee = 0 
(A71) 
and multiplying Eq. (A67) by A gives 
aS cae 
A aan Ans = 0 (A72) 


On adding Egs. (A71) and (A72) to Eq. (A644), and introdugima7 


Iq _ cua 
. er hw oe 
and 
aA GA 
A Ss Js - A §S a6 = @) (A yas) 











Od JS JA ore JS aA aq 
0 Ee eee es 9 as 
2 9S aA _ yl oq oA oq 
+ A Bre eee as ae 9) As Sire oS ee el ies 
2A aA y-l 34 ab, 
7) fas) See qASIa— + cos 6 ee = 0 
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Rearranging and introducing Eq. (A69) gives 


2Q wy = weet I |e 
ie SSS Sidaac air ya a US 21 Rymasei aE) | 
_ S qaA S (22 +cos 6 oo) (A75) 


Epeeeubtracting Eqs. (A71) and (A72) from Eq. (A64), again using 


Peet A?S) and” Eq. (A744), and introducing Eq. (A700), one also 


obtains 
3R 4) §R 2 yz 27/9. 2 
ee a ee AIS Tyg ty AN 
+ aq as(2? + cos 0 So) (A76) 


Megecther with Eqs. (A65) through (A67), Eqs. (A75) and (A76) 
form the system of equations which describe the isentropic flow 
of an inviscid perfect gas under unsteady, compressible condi- 
tions, and which may be solved as a system of ordinary differ- 
ential equations along characteristic trajectories in the 


space-time domain. 


ie) seOCK JUMP EQUATION 

The analytical expression for the change in the extended 
Riemann variable, Q, across a stationary shock is derived 
below. 

Consider a shock moving with a velocity Ve Tol xe) ela wipdlyelatys| 


with velocity dp as depleted etmieicemrere With the conditions 
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S 
a. ame > dp en mes it U; - doa 
A B A B 


Figure Al2 Figure Al3 


to the left of the shock being denoted by A subscripts and 
conditions to the right by B subscripts. If a velocity of ae 
is imposed on the entire system, the situation becomes one of 
a stationary shock, as depicted in Fig. Al3. In both cases, 
the high pressure side is to the left (A side). Since all 
velocities are defined positive to the right, a positive value 


for the relative mach number is defined as 





W (q, -V_) 
wWo= - me eed + s (A77) 
B B 
The extended Riemann variable, Q, 1S defined as 
Oieg= ol ChawAGS (A78) 
where 
es 2 al p 
= = ‘ g cae A 
° yl 7 Yey-ty Py ae 


p 


Tf all velocities are non-dimensionalized by A] and prescu 


and densities are non-dimensionalized by their values at B, 
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tat Aa "A, : Gp ts "B, 
AR AR 

SA “A Se oS 

ne An A B 
eAem “Al 2 ee ay 

_ Ap Y-L ¥(y-1) “pp op 

i Lp, SB (By "y 

¥oL. 9 yee) Pp Pp 
fa Se Ay 2 (Ap pen PA(PBYY, 

AG Ar i An ¥(y-1) Pp PA 

(A80) 


The ratios of pressure, density and sound speed across a 


normal shock are well known functions of yY and the Mach number 


meh in this case is W. 


An expression 


conservation. 


They are 
aE a esi (A81) 
(yt) we (A82) 
Coote ee 
eee ey ee (A83) 


fs _ ; 
Les (qn ee) as is obtained using mass 
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AA lavaea 
MA PBL yeh +2 
“B Al (yt+l) w 


Subtracting one from both sides, 


Ya TYR _ tyoL)w +2 = (y#l) Ww 
“B (y+1) we 
Ue =. 2 
ee oie 2 OL (A84) 
B (oy eet 
Substituting Ba.) vay? )) few Up, in Eq. (A84) gives 
Vata See ene 
7. <i 
B Cy) aw 
u -u, _ aa =a) nee 
A (vy 
B 
Since Uy = oe and Up = dp -V. 
Un Ups Sone cue one an © nem, 


Thus Eq. (A85) can be written as 


ae 


Sorts 7 2(w* <1) 


Ae ~  (yt+l)w (A86) 


Substituting Eqs. (A86) and (A81) through (A83) into Eq. (A80) 


gives 
= ; M/Z 
A B 2(w -l1) 2 (2 eee 
Ag (y+l1) W i Ta Gaal Oo L) (1 sea eae 14 
2 - a idee). 2 ieee ha ie 2 
iL "(Gee 1) w ) (aw 1d) } 


(A87) 
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(y+1l)w 
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APPENDIX B 


USING BULER=1) CN Inter o yi CiSeSs sie 


A. MEMORY REQUIREMENTS 

Storage should be defined as one mega-byte to ensure ade- 
quate memory for DISSPLA requirements. Also ensure that room 
exists on the user's disk for the listing file if that 


output option is selected. 


B. TERMINAL REQUIREMENTS 

The type of terminal required to run EULER-1l depends on 
the output option selected. When using the tabular listing 
output, any terminal tied into the VM/CMS system mav be 
used. Graphical output may also be created uSing any terminal 
and the graph stored on the user's disk ina metafile for later 
viewing at a graphics terminal or for printing. To have the 
graph stored on the disk, use the "COMPRS" command on line 
223 of the vrogram and comment out the “*TEK618" conmand an 
line 224. When graphical output is selected and it is desired 
to view the graph at the terminal at execution time, an IBM 
3277-Tek618 dual screen terminal must be used. To use this 
option, use the "TEK618" command on line 224 and comment out 
the "COMPRS" command on line 223. When uSing this option, a 
low quality hard copy of the TEK618 display may be obtained 


using the TEK4631 hard copy unit attached. 


wz 


C. PROCEDURE AND COMMANDS 

After loging on to VM/CMS, use the command "DEFINE 
STORAGE 1M" to increase the virtual memory as recommended 
above. This will remove you from the CMS environment. Re- 
enter CMS with the command "I CMS." 

Edit the user input area of EULER-1 FORTRAN as desired 
for the particular problem being run. When graphical output 
is selected, comment out either the "COMPRS" or the "TEK618" 
command as desired. 

Senpute the program uSing the command “FORTVS EULER-1." 
After compilation, an EULER-1l TEXT and an EULER-1 LISTING 
file will reside on the user's disk. At execution, if a 
tabular output has been selected, it will be sent to the 
user's disk as "FILE FTO9F001." To give this file a particu- 
lar name, say "EULER-1l LISTING," use the file definition 


command 
"er T 9 DISK EULER-1 LISTING A ( PERM" 


at compile time. An alternate and convenient means of com- 
piling the program is to set up an exec file on the user's 


disk by the name "EULER-1l EXEC" which contains the commands 


Bete Dt BULER-_ BISTENG A ( PERM 


POR vs EUEER-1 


Then to compile the program and define the output file, use 
Pac ecommand =EULER-1.” 
eer wecemprtatlon, to execute the program, use the 


command "DISSPLA EULFR-1." The message 


3) 3 


»-- YOUR FORTRAN PROGRAM IS NOW BEING LOADED... 


---BXECUTION WILEB SOON DROLLOw a 
should appear, followed shortly by 
72S AeeU TPO monG UN Sees 


If tabular output was selected, proper termination will 
result in a ready message. If graphical output was selected 
using the TEK618 option, the plot will appear on the TEK618 


terminal, the program will pause, and a 
Uo «sPLess “ENTER = te Ween nue § 


message will be displayed on the 3277 terminal. At this 
point the hard copy must be made, if desired, using the 
TEK4631 hard copy unit. After pressing the ENTER key on the 
3277 terminal, the plot will be erased and the program will 


terminate. Proper termination will be indicated by the message 
"END OF @DISSPEAY Voz HHH? VECTORS IN 1 PLOT 23a 


and a ready message. 
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